The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48_1/new2q.wic

Internet Message Format | 1995-03-23 | 5KB

From comp.sys.handhelds Thu Apr  4 13:04:21 1991
Path: uncw!ecsvax!ecsgate!mcnc!uvaarpa!haven!uflorida!caen!sdd.hp.com!hp-pcd!hpcvra.cv.hp.com!billw
From: [e-mail]. (William C Wickes)
Newsgroups: comp.sys.handhelds
Subject: ->Q or Not ->Q
Message-ID: <[e-mail].>
Date: 3 Apr 91 01:07:32 GMT
Organization: Hewlett-Packard Co., Corvallis, OR, USA
Lines: 108

From the HP 48 Design Team:


               ->Q OR NOT ->Q, THAT IS THE QUESTION

Thanks to Joseph Horn, the HP28/48 community now has additional fast
rationalization functionality for their machines (see "HP 48 Improved
->Q".)  We congratulate him for this elegant solution!  There are a
number of interesting points regarding DEC2FRACs relation to ->Q
which deserve some mention.  Perhaps these will stimulate discussion
and further contributions.

DEC2FRAC is an addition and not a replacement for ->Q because its aim
is rather different. DEC2FRAC's aim is to produce the fraction with
the _smallest error_ and a given denominator size. ->Q's aim is to
produce the fraction with the _smallest denominator_ and a given error
size. This is best illustrated with the golden mean, '(1+\sqr(5))/2'.
Given this and 1E11, DEC2FRAC returns '139583862445/86267571272' while
->Q returns (in STD mode) '514229/317811'. When evaluated, these
return _exactly_ the same floating point number. The idea of ->Q is
to "round to simpler" rather than "round to closer." This is useful
where you have what you believe to be a "simple" fraction contaminated
with roundoff or other noise.

The HP48 manuals did little to clarify this issue, especially since we
wanted to make no specific claim with only a conjecture of correctness. In
this regard, posting a proof that the algorithm fills in all possible
fractions would be a great boon to the community.

While for many "interesting" numbers the continued fraction sequence as
generated by DUP IP SWAP FP INV ... is exact, even though the floating
point representation is not, in other cases, you get "noise" unexpectedly
early. Does this have any significance regarding the "safety" of generating
the numerator of the fraction from the denominator by division?

While ->Q keeps around the entire continued fraction sequence entailing
some speed penalty, there may be other statistics about the sequence
(such as periodicity) which can be of interest in "deciphering" a
floating point number. Suggestions along these lines are most welcome.

Thanks again to J. Horn for a stimulating and useful post.


Included below is a program, NEW2Q which follows Horn's algorithm, but
uses exit conditions like those of ->Q.

@ NEW2Q: Version of ->Q based on J.K.Horn's Algorithm
@
@ Input:
@
@ 2: Decimal Number to be converted to a fraction
@ 1: Number of decimal places of zeros required in the error.
@
@ Output:
@
@ 1: 'Numerator/Denominator'
@
@ Example:
@
@ What's the simplest fraction which agrees with sqrt(3) to 4 decimal places?
@   3 \sqr 4 NEQ2Q returns '97/56'
@   '97/56-\sqr(3)' EVAL returns .00009294957
@                                 ^^^^  note 4 zeros.
@
@ Checksum and bytes:
@   #3992d
@   620.5

%%HP: T(3)A(R)F(.);
\<< DUP2
  IF 1 > SWAP FP AND
  THEN OVER XPON 1 -                        @ calculate the input exponent.
    \<< \-> X 'IFTE(X==0,-500,XPON(X))' \>> @ define a 0-tolerant XPON.
  \-> f c x expo
    \<< 0 1 1 f DUP IP SWAP FP              @ set recursion initial cond.s.
      WHILE
       OVER 5 PICK / f - ABS expo EVAL      @ calculate expon. of error
       x SWAP - c <                         @ and compare with input.
       OVER AND                             @ if not close enough and
                                            @ the remainder's non-zero
      REPEAT
       INV DUP FP SWAP IP                   @ then calculate next iterate.
       \-> B0 B1 A0 A1 R B
        \<< B1 'B*B1+B0' EVAL
            A1 'B*A1+A0' EVAL
            R
        \>>
      END
      DROP SWAP DROP SWAP
      DUP 4 ROLL - DUP f * 0 RND            @ calc. "missing" den. and num.
      \-> N D D0 N0
      \<<
        IF 'x-expo(ABS(f-N0/D0))<c'         @ if "missing" frac. is not
        THEN N D                            @ good enough, use original.
        ELSE N0 D0
          IF 'N0\=/N'                       @ If it is really new,
          THEN 200 .2 BEEP                  @ then beep.
          END
        END
      \>>
    \>>                                     @ We're done, now clean up.
    IF DUP ABS 1 >
    THEN # 352318d SYSEVAL
    ELSE DROP
    END
  ELSE DROP
  END
\>>